{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Шахова Нонна – @nonna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python R Collaboration\n",
    "### life hacking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Не знаю как вы, но я постоянно встречаю сравнение R и Python в задачах работы с данными и машинного обучения. Я бы даже сказала, их противопоставление и выяснение, какой же лучше? Вы можете найти массу исследований, рассматривающих практически все аспекты работы с тем и другим, но вряд ли где-то есть ответ: \"Используй один из двух и все будет OK вообще всегда!\" Скорее вы встретите что-то вроде: \"Для этой конкретной задачи лучше использовать R, а для той, конечно, Python\" - и тут не только в самой задаче дело. Очевидно, что все зависит от конечной цели. Нужно глубокое научное исследование? Бери R, там есть все для этого. Хочешь, чтоб работало быстро и в продакшене? Возьми лучше Python, потом спасибо скажешь! И так до бесконечности...\n",
    "\n",
    "<img src=\"../../img/RPy.jpg\">\n",
    "\n",
    "Предположим, у нас все написано на R либо Python: все хорошо работает и проблем не приносит. Но вдруг возникла гениальная мысль, которую хочется срочно проверить, а подходящего пакета нет?!\n",
    "\n",
    "Для этого существуют пакеты несколько способов взять все лучшее из Python и R и собрать своего Франкинштейна. Данный туториал будет об одном из таких способов решения самых тривиальных задач, с которыми вы сталкивались или обязательно повстречаетесь в будущем:\n",
    "\n",
    "\n",
    ">- Как определить число кластеров и не ошибиться?\n",
    ">- Как представить в 2D на данные с большим количеством признаков и не потерять суть?\n",
    "\n",
    "Прежде чем начать работу:\n",
    "\n",
    "1 способ: Если вы не используете докер:\n",
    "    - поставьте R под вашу систему (OS X, Win, Ubu): https://www.r-project.org \n",
    "    - поставьте rpy2: conda install rpy2\n",
    "    - остальные R пакеты потренируемся ставить из jupyter notebook\n",
    "    \n",
    "   \n",
    "2 способ: Если вы используете докер-контейнер:\n",
    "    - придется поставить R в контейнер и собрать образ снова (т.к. лозунг этого туториала \"по-быстренькому\", если вы новичок, который ранее не собирал контейнеры, следуйте первому варианту)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import rpy2.interactive as r\n",
    "import rpy2.interactive.packages\n",
    "from rpy2.robjects import pandas2ri\n",
    "from sklearn import datasets\n",
    "\n",
    "pandas2ri.activate()\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "rlib = r.packages.packages\n",
    "r.packages.importr(\"utils\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Лафхак 1: подбираем число кластеров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предположим, что мы исследуем вопрос об оптимальном количестве кластеров для какого-либо датасета. Общеизвестное правило локтя покажет нам такую картину:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../../img/Noclast.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "По такому графику сложно сказать уверенно, какое же количество кластеров будет лучшим. Кто-то скажет 2, кто-то 4 - крайне субъективная оценка... И тут на помощь приходит мудрость толпы - устроим голосование!\n",
    "\n",
    "Для этого используем 30 показателей оптимального количества кластеров R пакета NbClust: https://cran.r-project.org/web/packages/NbClust/NbClust.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Для примера возьмем всем хорошо знакомый датасет Iris\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "df = pd.DataFrame(data = iris.data, columns = iris.feature_names)\n",
    "\n",
    "features = df.columns\n",
    "df['target'] = iris.target\n",
    "import seaborn as sns\n",
    "\n",
    "%pylab inline\n",
    "\n",
    "sns.pairplot(df[features]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Визуально можно выделить два класса, но мы подозреваем, что их больше, так сколько же?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# преобразуем его в R-пригодный вид\n",
    "r_df = pandas2ri.py2ri(df[features])\n",
    "\n",
    "# и отправим в R\n",
    "%Rpush r_df\n",
    "\n",
    "# теперь пришло время поставить в R пакет NbClust, сделать это можно так\n",
    "rlib.utils.install_packages(\"NbClust\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# а теперь начинается магия!\n",
    "\n",
    "# подключаем библиотеку\n",
    "%R library(NbClust)\n",
    "\n",
    "# сид для воспроизводимости\n",
    "%R set.seed(1234)\n",
    "\n",
    "# и указываем диапазон от 2 до 10 кластеров\n",
    "%R Nnc <- NbClust(r_df, distance = \"euclidean\", min.nc = 2, max.nc = 10, method = \"complete\", index = \"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Цитируя Conclusion: ...the best number of clusters is  3! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# интересно посмотреть как распределились голоса\n",
    "%R table(Nnc$Best.n[1,])\n",
    "%R barplot(table(Nnc$Best.n[1,]), xlab=\"Number of Clusters\", ylab=\"Number of Criteria\", main=\"Number of Clusters Chosen by 30 Criteria\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выглядит очень убедительно, давайте посмотрим на данные в привычном виде - нет ли явных противоречий?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(df, hue = 'target');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Очевидно, что данные вполне хорошо разделимы на три кластера, как и выяснилось с помощью NbClust.\n",
    "\n",
    "На самом деле, из этого пакета можно вытащить намного больше полезной информации, особенно это может быть полезно в Research работах."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Лафхак 2: организуем самоорганизующуюся карту Кохонена за 5 мин"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В двух словах, самоорганизующиеся карты Кохонена славятся тем, что делают 2D представление из многомерного(и даже очень) представления с минимумом потерь смысла, а еще они самообучаются. R-пакет kohonen обладает массой полезных инструментов, например, heatmaps и разрезы по каждому признаку. А также это очень красиво!\n",
    "\n",
    "https://cran.r-project.org/web/packages/kohonen/kohonen.pdf\n",
    "\n",
    "В данном примере возьмем датасет побольше: Boston data (недвижимость). В данном случае, мы не решаем задачу кластеризации, но хотим использовать аггрегирование данных как дополнительную feature. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = datasets.load_boston()\n",
    "b_df = pd.DataFrame(data = boston.data, columns = boston.feature_names)\n",
    "\n",
    "br_df = pandas2ri.py2ri(b_df)\n",
    "sns.pairplot(b_df);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# поставим в R пакет kohonen\n",
    "rlib.utils.install_packages(\"kohonen\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%Rpush br_df\n",
    "\n",
    "%R library(kohonen)\n",
    "%R set.seed(1234)\n",
    "\n",
    "# данные нужно отдать в виде матрицы\n",
    "%R b_matrix <- as.matrix(br_df)\n",
    "%R som_grid <- somgrid(xdim = 20, ydim=20, topo=\"hexagonal\")\n",
    "\n",
    "# и организовать карту\n",
    "%R bos.som <- som(b_matrix, som_grid, rlen=100, alpha=c(0.05,0.01));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# теперь определим количество кластеров уже известным нам трюком\n",
    "%R codes <- getCodes(bos.som)\n",
    "%R b_Nnc <- NbClust(codes, min.nc=2, max.nc=15, method=\"kmeans\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# и посмотрим на некоторые визуализации\n",
    "%R plot(bos.som, type=\"dist.neighbours\")\n",
    "%R plot(bos.som, type=\"changes\",main=\"Training Progress\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# в конечном итоге, вот так выглядит в двумерном пространстве весь датасет\n",
    "%R som_cluster <- cutree(hclust(dist(codes)), 3)\n",
    "%R pretty_palette <- c('#7DBD00', '#62C2CC', '#FF5B00', '#f1bc41')\n",
    "%R plot(bos.som, type='mapping', bgcol = pretty_palette[som_cluster], main = \"Clusters\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если немного поколдовать над визуальной составляющей, то можно получать более приятные изображения, например:\n",
    "<img src=\"../../img/NCl.png\">\n",
    "Подводя итог, хочется отметить, что самоорганизующиеся карты удобный, понятный, но незаслуженно обойденный питоном инструмент. Я пробовала, наверное, все реализации этого пакета для питона, но они существенно уступают R. Конечно, карты Кохонена будут использоваться не каждый день, но, согласитесь, это очень удобно: запихать туда весь датафрейм, оценить, насколько четкая или нет идет кластеризация и принять решение, имеет ли смысл вводить такой признак или сэкономить время для другой полезной работы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Что осталось за кадром"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Эти два крошечных примера не могут показать всех преимуществ использования R и Python, а они все больше становятся взаимопомогающими, а не противоборствующими. Можно привести пример в той же ARIMA, которая настолько прекрасна в R, что переписывать ее на Python нет большого смысла. Ведь к ней также можно обращаться, особенно для подбора p, d, q, P, D, Q:\n",
    "> %R model_lam = forecast::auto.arima(ts_a, stepwise = TRUE, lambda = ll)\n",
    "\n",
    "https://cran.r-project.org/web/packages/forecast/forecast.pdf\n",
    "\n",
    "А также выделить тренд из STL (который часто бывает полезной feature)\n",
    "> %R trendSTL = stl$time.series\n",
    "\n",
    "https://cran.r-project.org/web/packages/stlplus/stlplus.pdf\n",
    "\n",
    "Все это очень удобно для первого приближения к желаемой модели, т.к. где-то в идеальном мире мы не используем пакеты из коробки, а детально разбираемся в том, как каждый из них работает ;)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
